%% LDPC + M-QAM + AWGN/Rayleigh 彩色图像传输仿真 (优化版)
% 说明：
%   - 读取彩色 (RGB) 图像，将其转为比特流（24 bit）；
%   - 使用 LDPC 1/2 码进行编码；
%   - 对编码后比特流进行 M-QAM 调制 (支持 16, 32, 64)；
%   - 通过 AWGN 或 Rayleigh 信道传输，可设置任意 SNR 列表；
%   - 使用向量化/并行计算减少循环，提升运行速度；
%   - 解调、解码、重构彩色图像，并绘制 BER 曲线，以及展示恢复图像。
% 
% 要求：
%   - 需安装 Communications Toolbox；
%   - 如可使用 Parallel Computing Toolbox，可启用 parfor。

clear; clc;

%% ==== 可配置参数 ====
ldpcRate      = 1/2;                    % LDPC 码率：1/2
qamOrders     = [16, 64];               % 常用 QAM 阶数，可改为 [16 32 64]
channelTypes  = {'AWGN', 'Rayleigh'};    % 可选 'AWGN' 或 'Rayleigh'
snrVec_dB     = 0:2:16;                 % SNR 列表 (dB)
maxLDPCIter   = 50;                     % LDPC 最大迭代次数
useParallel   = true;                   % true: 尝试并行计算（需 Parallel Toolbox）
showImage     = true;                   % 是否展示恢复图像（仅第一个 SNR）

%% ==== 读取并预处理彩色图像 ====
imgRGB = imread('peppers.png');         % 可替换成其他彩色图像
if size(imgRGB,3) ~= 3
    error('输入图像必须为彩色 RGB 图像');
end
[Himg, Wimg, C] = size(imgRGB);
numPixels    = Himg * Wimg;
fprintf('原始彩色图像：%d × %d × %d = %d 像素 (每像素 3 通道)\n', Himg, Wimg, C, numPixels);
% 转为 (numPixels × 3) 的矩阵
pixelsRGB    = reshape(imgRGB, [], 3);       % 每行 [R G B] (0~255)
% 对 R/G/B 分别拆为 8-bit
bitsR        = de2bi(pixelsRGB(:,1), 8, 'left-msb');  % (numPixels×8)
bitsG        = de2bi(pixelsRGB(:,2), 8, 'left-msb');  % (numPixels×8)
bitsB        = de2bi(pixelsRGB(:,3), 8, 'left-msb');  % (numPixels×8)
% 拼接成 (numPixels × 24) 的比特矩阵，顺序为 [R7..R0, G7..G0, B7..B0]
bitsRGB      = [bitsR, bitsG, bitsB];       % (numPixels×24)
% 转成一维比特流（按行展开）
txBits       = reshape(bitsRGB.', [], 1);    % (numPixels*24×1)
lenInfoBits  = numel(txBits);

%% ==== LDPC 校验矩阵加载 (IEEE 802.11n rate-1/2) ====
H = ldpcParityCheckMatrix('802.11n','Rate1/2');
nLDPC = size(H,2);
kLDPC = nLDPC * ldpcRate;
fprintf('使用 LDPC (n=%d, k=%d) 码，码率=%.2f\n', nLDPC, kLDPC, ldpcRate);

ldpcEnc = comm.LDPCEncoder(H);
ldpcDec = comm.LDPCDecoder(H, ...
    'IterationTerminationCondition','MaxLLR', ...
    'MaximumIterationCount', maxLDPCIter);

%% ==== 信息块分割与 LDPC 编码（预处理，只执行一次） ====
numBlocks   = ceil(lenInfoBits / kLDPC);
padBits     = numBlocks * kLDPC - lenInfoBits;
txBitsPad   = [txBits; false(padBits,1)];  % 零填充到 kLDPC
% 预分配编码后比特数组
txBitsLDPC  = false(numBlocks * nLDPC, 1);
for blk = 1:numBlocks
    i1 = (blk-1)*kLDPC + 1;
    i2 = blk*kLDPC;
    infoBlk = txBitsPad(i1:i2);
    txBitsLDPC((blk-1)*nLDPC + (1:nLDPC)) = ldpcEnc(infoBlk);
end
lenTxLDPC = numel(txBitsLDPC);
disp(['LDPC 编码后比特长度：', num2str(lenTxLDPC)]);

%% ==== 对每种 QAM 阶数逐一仿真 ====
for qIdx = 1:numel(qamOrders)
    M = qamOrders(qIdx);
    bitsPerSym = log2(M);
    fprintf('\n===== %d-QAM 仿真 =====\n', M);
    % QAM 调制器/解调器 对象（预创建，提高效率）
    qamMod = comm.RectangularQAMModulator( ...
        'ModulationOrder', M, ...
        'BitInput', true, ...
        'NormalizationMethod', 'Average power');
    qamDemod = comm.RectangularQAMDemodulator( ...
        'ModulationOrder', M, ...
        'BitOutput', true, ...
        'DecisionMethod', 'Approximate log-likelihood ratio', ...
        'NormalizationMethod', 'Average power');

    % QAM 前填零：保证比特数是 bitsPerSym 的整数倍
    padQAM = mod(-lenTxLDPC, bitsPerSym);
    if padQAM < 0, padQAM = padQAM + bitsPerSym; end
    bitsForQAM = [txBitsLDPC; false(padQAM,1)];
    symNum     = numel(bitsForQAM) / bitsPerSym;
    txSyms      = qamMod(bitsForQAM);  % 一次性调制

    for cIdx = 1:numel(channelTypes)
        chanType = channelTypes{cIdx};
        fprintf('\n-- 信道类型：%s --\n', chanType);
        % 创建 Rayleigh 信道对象（如需）
        if strcmpi(chanType,'rayleigh')
            rayChan = comm.RayleighChannel( ...
                'SampleRate',1, ...
                'PathDelays',0, ...
                'AveragePathGains',0, ...
                'NormalizePathGains',true, ...
                'MaximumDopplerShift',0);
        end
        % 预分配 BER 向量
        berVec = zeros(size(snrVec_dB));
        % 选择循环方式：parfor 加速
        usePar = useParallel && license('test','Distrib_Computing_Toolbox');
        if usePar
            parfor sIdx = 1:length(snrVec_dB)
                berVec(sIdx) = simulateSNR(snrVec_dB(sIdx), chanType, txSyms, qamDemod, ldpcDec, numBlocks, nLDPC, kLDPC, lenInfoBits, rayChan, txBits);
            end
        else
            for sIdx = 1:length(snrVec_dB)
                berVec(sIdx) = simulateSNR(snrVec_dB(sIdx), chanType, txSyms, qamDemod, ldpcDec, numBlocks, nLDPC, kLDPC, lenInfoBits, rayChan, txBits);
            end
        end

        % 绘制 BER 曲线
        figure;
        semilogy(snrVec_dB, berVec, '-s','LineWidth',1.2);
        grid on;
        xlabel('SNR (dB)'); ylabel('BER');
        title(sprintf('LDPC+%d-QAM over %s 信道 BER', M, chanType));

        % 恢复并展示第一个 SNR 对应的彩色图像
        if showImage
            SNR0 = snrVec_dB(1);
            recImg = recoverColorImageAtSNR(SNR0, chanType, txSyms, qamDemod, ldpcDec, numBlocks, nLDPC, kLDPC, lenInfoBits, rayChan, Himg, Wimg, txBits);
            figure;
            subplot(1,2,1); imshow(imgRGB); title('原始彩色图像');
            subplot(1,2,2); imshow(recImg); title(sprintf('%d-QAM, %s, SNR=%d dB', M, chanType, SNR0));
        end
    end
end

fprintf('仿真完成。\n');

%% ==== 子函数：在给定 SNR 下仿真并返回 BER ====
function ber_val = simulateSNR(SNR, chanType, txSyms, qamDemod, ldpcDec, numBlocks, nLDPC, kLDPC, lenInfoBits, rayChan, txBits)
    if strcmpi(chanType,'awgn')
        rxSym = awgn(txSyms, SNR, 'measured');
    else % Rayleigh + AWGN
        faded = rayChan(txSyms);
        rxSym = awgn(faded, SNR, 'measured');
        reset(rayChan);
    end
    rxLLR = qamDemod(rxSym);
    % 预分配解码后的比特存储
    decBits = false(numBlocks * kLDPC, 1);
    % 按块 LDPC 解码
    for blk = 1:numBlocks
        idxL = (blk-1)*nLDPC + (1:nLDPC);
        llrBlk = rxLLR(idxL);
        decBits((blk-1)*kLDPC + (1:kLDPC)) = ldpcDec(llrBlk);
    end
    decBits = decBits(1:lenInfoBits);
    ber_val = sum(decBits ~= txBits) / lenInfoBits;
end

%% ==== 子函数：恢复彩色图像在给定 SNR ==== 
function recImg = recoverColorImageAtSNR(SNR, chanType, txSyms, qamDemod, ldpcDec, numBlocks, nLDPC, kLDPC, lenInfoBits, rayChan, Himg, Wimg, txBits)
    if strcmpi(chanType,'awgn')
        rxSym = awgn(txSyms, SNR, 'measured');
    else
        faded = rayChan(txSyms);
        rxSym = awgn(faded, SNR, 'measured');
        reset(rayChan);
    end
    rxLLR = qamDemod(rxSym);
    decBits = false(numBlocks * kLDPC, 1);
    for blk = 1:numBlocks
        idxL = (blk-1)*nLDPC + (1:nLDPC);
        decBits((blk-1)*kLDPC + (1:kLDPC)) = ldpcDec(rxLLR(idxL));
    end
    decBits = decBits(1:lenInfoBits);
    % 还原比特到像素矩阵 (numPixels × 24)
    bitsMat = reshape(decBits, 24, []).';
    % 分离 R, G, B 通道各 8-bit
    R_bits = bitsMat(:,1:8);
    G_bits = bitsMat(:,9:16);
    B_bits = bitsMat(:,17:24);
    R_vals = bi2de(R_bits, 'left-msb');
    G_vals = bi2de(G_bits, 'left-msb');
    B_vals = bi2de(B_bits, 'left-msb');
    % 合并并 reshape 回 Himg×Wimg×3
    pixelsRec = uint8([R_vals, G_vals, B_vals]);
    recImg = reshape(pixelsRec, Himg, Wimg, 3);
end
